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Abstract: Recent studies on the Lorentzian version of the type IIB matrix model show 
that (3+l)D expanding universe emerges dynamically from (9+l)D space-time predicted 
by superstring theory. Here we study a bosonic matrix model obtained by omitting the 
fermionic matrices. With the adopted simplification and the usage of a large-scale par¬ 
allel computer, we are able to perform Monte Carlo calculations with matrix size up to 
N = 512, which is twenty times larger than that used previously for the studies of the 
original model. When the matrix size is larger than some critical value W ~ 110, we 
find that (3-1-1 )D expanding universe emerges dynamically with a clear large-iV scaling 
property. Furthermore, at sufficiently late times, we observe a power-law behavior of 
the spatial extent with respect to time t, which is reminiscent of the expanding behavior 
of the Friedmann-Robertson-Walker universe in the radiation dominated era. We discuss 
possible implications of this result on the original model including fermionic matrices. 
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1. Introduction 

Understanding how our universe began is a fascinating subject in theoretical physics. This 
is, of course, not so easy since one has to deal with quantum gravity near cosmic singular¬ 
ity, which appears in general relativity. As a promising way to describe quantum gravity, 
string theory has been studied for many years. However, as far as one applies pertur¬ 
bative string theory around a time dependent background, cosmic singularity cannot be 
resolved in general [1-4]. It is therefore necessary to study string theory in a nonperturba- 
tive background-independent fashion. As nonperturbative formulations of superstring/M 
theory, matrix models that can be obtained formally by dimensionally reducing lOD J\f = 1 
super Yang-Mills theory to d = 0 [5], d = 1 [6] and d = 2 [7] were proposed. As a closely 
related direction, ref. [8,9] proposes a conformal field theory, which is holographically dual 
to inflationary models. 

The type IIB matrix model [5] is one of these proposals corresponding to the d = 0 case 
above. A peculiar feature of this model is that both space and time emerge dynamically 
from the matrix degrees of freedom. In this context, the idea of emergent gravity has been 
pursued [10-17] in gauge theories on non-commutative space that appear from the type IIB 
matrix model for a particular class of backgrounds [18-21]. There are also other proposals 
for the description of curved space-time in matrix models [22-24]. 

Until recently, the type IIB matrix model was studied after the “Wick rotation” [25-38] 
since the partition function of the Euclidean model obtained in this way is shown to be 
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finite [39,40]. However, the Euclidean model is clearly not suitable for studying the real¬ 
time dynamics since the time coordinate is treated as purely imaginary due to the “Wick 
rotation”. Moreover, it is known that the Wick rotation is more subtle in quantum gravity 
than in quantum field theories at the nonperturbative level (See, for instance, refs. [41,42]). 
In fact, according to a recent study of the Euclidean model using the Gaussian expansion 
method, the emergent space-time seems to have only three dimensions [43]. 

On the other hand, the Lorentzian version of the type IIB matrix model has been 
studied for the first time in ref. [44]. Unlike the Euclidean version, one has to introduce 
infrared cutoffs in the temporal and spatial directions to make the partition function finite. 
Despite this subtlety, the time-evolution of the “universe” was extracted from the matrix 
configurations generated by Monte Carlo simulation, and a large-scaling behavior was 
observed with the matrix size N < 16. Quite remarkably, the SO(9) rotational symmetry 
of the 9D space turned out to be spontaneously broken down to SO(3) at some critical 
time, after which only three out of the nine spatial directions start to expand. 

In order to see what happens at later times in this model, one needs to increase the 
matrix size, which makes the Monte Carlo simulation more and more time-consuming. 
As an alternative approach, one can use the classical approximation [45-49] to investigate 
possible behaviors at late times since each term in the action becomes large as the expansion 
of the “universe” proceeds. A general prescription to construct solutions to the classical 
equations of motion was given in ref. [46]. One can actually construct classical solutions 
corresponding to an expanding (3-|-l)D universe, which naturally solve the cosmological 
constant problem [46]. As a closely related progress, it was found that matrix configurations 
with intersecting fuzzy spheres in the extra dimensions can accommodate the standard 
model fermions [50-56]. 

In fact, it is known that the classical equations of motion of the matrix model have 
infinitely many solutions [46]. Therefore, in order to determine which classical solution is 
actually realized at late times, we need to study the time-evolution of the “universe” at 
least for a sufficiently long time by performing Monte Carlo simulation. To that end, we 
previously studied a simplified model that describes the early time behaviors of the original 
model [57]. With the matrix size N < 64, we observed a clear exponentially expanding 
behavior, which is reminiscent of the inflation. Monte Carlo studies of the original model 
with A^ = 24 [58] yielded results consistent with this observation. 

In this paper we study a bosonic model, which can be obtained by simply omitting 
the fermionic matrices. This simplification and the usage of a large-scale parallel computer 
enable us to perform Monte Carlo simulation with the matrix size up to A^ = 512. Unlike 
the original model, the eigenvalue distribution of the temporal matrix turns out to have 
a finite extent without introducing a cutoff in the temporal direction. The extent of the 
eigenvalue distribution is independent of N for N < Nc — 110, but it increases with N 
for N > Nc- We find that the properties of the model changes drastically at the critical 
N = Nc- For N < Nc, the dominant matrix configurations do not allow extraction of a 
well-dehned time-evolution. For N > Nc, on the other hand, we can extract a meaningful 
time-evolution, which shows that the SO(9) rotational symmetry is broken spontaneously 
down to SO(3) symmetry at some point in time similarly to the original model. The large- 
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N scaling behavior is clearly observed. The expanding behavior of the spatial extent can be 
fitted by an exponential function only for a rather short period, and after that it becomes 
a power-law with respect to time t. The latter behavior coincides with the expanding 
behavior of the Friedmann-Robertson-Walker universe in the radiation dominated era. 
From these results, we speculate that the exponential expansion of the space in the original 
model suggested in the previous works actually terminates at some point in time and turns 
into a power law similarly to the bosonic model. This would imply that the number of 
e-foldings is determined dynamically in the Lorentzian type IIB matrix model. 

The rest of this paper is organized as follows. In section we briefly review some 
important properties of the Lorentzian type IIB matrix model. In section ^ we define 
the bosonic model we study in this paper, and discuss the existence of the critical N, at 
which the properties of the model change drastically. In sections ^ and we discuss in 
detail the properties of the model below and above Nc, respectively. Section is devoted 
to a summary and discussions. In Appendix ^ we give the details of our simulation. In 
Appendix]^ we present our results for the (5-|-l)D version of the bosonic model. 


2. Brief review of the Lorentzian type IIB matrix model 

The action of the Lorentzian type IIB matrix model is given by [5] 


5 = 5b + 5f 

1 

4c?2 


5b = —Tr([A^,A,][A^A"]) , 




( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 


where the bosonic N x N matrices (/r = 0,..., 9) and the fermionic matrices To 
(« = !,..., 16) are both traceless and Hermitian. F^ are lOD gamma-matrices after the 
Weyl projection and C is the charge conjugation matrix. The “coupling constant” g is 
merely a scale parameter in this model since it can be absorbed by rescaling and T 
appropriately. The indices fj, and are contracted using the Lorentzian metric = 
diag (—1, !,...,!). The Euclidean version can be obtained by making the “Wick rotation” 
Aq = fAio, where Aio is supposed to be Hermitian. 

The partition function for the Lorentzian version is proposed in ref. [44] as 


Z = j dAd^ 


(2.4) 


with the action (13). The “i” in front of the action is motivated from the fact that the 
string world-sheet metric should also have a Lorentzian signature. By integrating out the 
fermionic matrices, we obtain the Pfaffian 


J d^ = PfA4 (A) , 


which is real unlike in the Euclidean case [38]. Note also that the bosonic action 
be written as 

i (Foif + rv } . 


(2.5) 


can 


( 2 . 6 ) 
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where we have introduced the Hermitian matrices = i Since the two terms 

in the last expression have opposite signs, 5b is not positive semi-definite, and it is not 
bounded from below. 

In order to make the partition function (2.4) finite, one needs to introduce infrared 
cutoffs in both the temporal and spatial directions, for instance, as [44] 


^Tv{Aof < K^Tr{Aif 

1 


N 


Tr (Aif < 


(2.7) 

( 2 . 8 ) 


It was found that the two parameters k and A can be removed in the large-limit in such a 
way that physical quantities scale [44]. Therefore the resulting theory can be characterized 
by a single scale parameter. 

In the present work, it will be important to understand the reason why we need to 


introduce the cutoff (2^) in the temporal direction. Note first that one can use the SU (N) 
symmetry of the model to bring the temporal matrix Aq into the diagonal form 


Aq = diag (ai,..., on) , where ai < • • • < on ■ 

By “fixing the gauge” in this way, we can rewrite the partition function (^(J) as 

r ^ r 

YldaaA{af / dAid^ 


Z = 


js 


a=l 


N 


A(a) = ]^ (tta - Ob) , 


(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


a>b 


where A(a) is the van der Monde determinant. The factor A(a)^ in ( 2.10] ) appears from the 
Fadeev-Popov procedure for the gauge fixing, and it acts as a repulsive potential between 
the eigenvalues at of Aq. Here we consider a situation in which the eigenvalues of Ao are 
well separated from each other. Then the action 5 = 5b -|- 5f can be expanded as 


5b — T'‘‘ ) 

‘S'f = -X-2 ('^'a)6a(aa — «&) (Cr^)c^y3 (^/3)afe 4- 


( 2 . 12 ) 

(2.13) 


omitting the sub leading terms for large [oa — afe]. Integrating out Ai at one loop neglecting 
the zero modes corresponding to diagonal elements, we obtain A(a)“^®. On the other hand, 
integrating out To at one loop similarly, we obtain A(a)^®. Thus we find that the potential 
between Oj is canceled exactly at the one-loop level. This is actually a consequence of 
supersymmetry [5] of the model (|2.4| ). Owing to this property, the eigenvalue distribution 
of Aq extends to infinity even for finite N if the cutoff ( p.7| ) were absent. 

In fact, after some manipulation and rescaling of A^, we can rewrite the partition 
function (2.4) as [44] (See Appendix A of ref. [57] for a rehned argument.) 


Z = / dAPfM 


(A) 6 (^Tr (T^.T/^")) 5 (^TV (Aif - l) 9 


K-lTr(Ao)2 
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(2.14) 


„ N a ^ ^ 

= / dAi A^{a) FfM (^) <5 f -TV 

a=l i=l ^ 

X ^ (;^Tr {Aif - 1^ 6» - ^Tr {A^f 

where 0 (x) is the Heaviside step function. This form allows us to performing Monte Carlo 
simulation of the Lorentzian model without the sign problem unlike the Euclidean model.^ 
It turns out that one can extract a time-evolution from configurations generated by 
simulating ( 2.141) . A crucial observation is that the spatial matrices Ai have a band-diagonal 
structure in the SU(A^) basis in which Aq has the diagonal form (|2.9| ). More precisely, there 
exists some integer n such that the elements of spatial matrices (Aj)^j, for \a — b\ > n are 
much smaller than those for |a — 5| < n. Based on this observation, we may naturally 
consider n x n matrices 

= , (2.15) 

as representing the state of the universe at time t, where I,J = 1,..., n and = 0,1,..., N— 


n. The time t in (2.15) is defined by 


1 ^ 

t = — > a^+i 
n 

i=\ 


(2.16) 


corresponding to the nxn matrices Aj. For example, we can define the extent of space at 
time t as 



(2.17) 


where the symbol tr represents a trace over the nxn block. We also define the “moment 
of inertia tensor” 

Tij (t) = ^tr (Ai (t) Aj (t)) , (2.18) 

which is a 9 X 9 real symmetric matrix. The eigenvalues of (t), which we denote by 
Aj (t) with the order 

Ai (t) > X 2 (t) > ■ ■ ■ > Xg (t) , (2A9) 


represent the spatial extent in each of the nine directions at time t. Note that the expec¬ 
tation values (Aj (t)) tend to be equal in the large-N limit if the SO(9) symmetry is not 
spontaneously broken. This is the case at early times of the time-evolution. After a critical 
time tc, however, we find that three largest eigenvalues (Aj (t)) (z = 1, 2, 3) become sig¬ 
nificantly larger than the others, which implies that the SO(9) symmetry is spontaneously 
broken down to SO(3). 

It would be interesting to study a long time-evolution of the model and see how the 
expansion of space proceeds. This requires very large matrices, which makes the simulation 
unfeasible. In the previous work [57], we studied a simplified model, in which the Pfaffian 

^Strictly speaking, the Pfaffian PfAl in ( ^.14| ) can change its sign, but it turned out that configurations 
with positive Pfaffian dominate at large N. 
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is replaced by the one-loop contribution A(a)^® mentioned above. This replacement is 
expected to be valid at early times, where the expansion of space has not proceeded much 
and the leading term in ( 2.13D is indeed dominant. According to the argument below 
( 2.131) , the potential between the eigenvalues of Aq is canceled at one loop and hence the 
cutoff ( p.Tj ) in the temporal direction is needed in this simplified model as well as in the 
original model. On the other hand, this simplihed model can be simulated with much less 
effort than the original model.^ In ref. [57] the (5+l)D version of the simplified model 
was studied with the matrix size N < 64, and the SO(5) symmetry was found to be 
broken spontaneously down to SO(3) at some point in time similarly to the original model. 
Moreover, the expanding behavior of the 3D space turned out to be exponential,^ and no 
tendencies of slowing down were observed within the scaling region. Analogous behaviors 
were also confirmed for the (9+l)D version of the simplified model. In the original model, 
on the other hand, the subleading term in the fermionic action ( 2.13| ) becomes important 
at late times as the expansion proceeds, and hence it can affect the expanding behavior. 


3. The bosonic Lorentzian type IIB matrix model 


In this paper we study a bosonic model in which the fermionic matrices are simply omitted. 
The partition function is given by 


Z = 


/ 




(3.1) 


In section ^ we reviewed an argument for the necessity of the temporal cutoff in the 
original model and the simplified model for early times. In the present case of the bosonic 
model o, the same argument implies that there is no need to introduce the temporal 
cutoff (^), and that we only need a cutoff (|2.8| ) in the spatial direction. Corresponding 
to (|t|), we can study the bosonic Lorentzian type IIB matrix model by simulating 


Z = I dA 6 (^TV 5 (^Tr (A,)' - l) (3.2) 

= , (3.3) 

a=l i=l ^ ^ ^ ^ 


which requires computational efforts comparable to the simplified model for early times 
reviewed in the previous section. We have used a large-scale parallel computer to simulate 
the model (|3.3| ) with the matrix size up to A = 512, which enables us to investigate a long 
time-evolution. See Appendix ^ for the details of the simulation. 

We measure the quantity (^Tr(Ao)^), which represents the extent of the eigenvalue 
distribution of Aq. As we mentioned above, this quantity turns out to be finite in the model 

^In order to make one trajectory in the Hybrid Monte Carlo algorithm, the original model requires 0{N^) 
arithmetic operations, whereas the simplified model requires only 0{N^) arithmetic operations. The reason 
for this is that the number of iterations required for the convergence of the conjugate gradient method used 
to implement the effects of fermions grows as 0{N^). 

®This behavior is also confirmed with smaller matrix size N < 32 with the aid of a renormalization group 
method developed in the same paper [57]. 
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Figure 1: (Left) The extent (^Tr(Ao)^) of the eigenvalue distribution of Aq is plotted against 
N. (Right) The expectation values \i {t) of the nine eigenvalues of T^- (t) at t = tpeak are plotted 
against N. For N < Nc — 112, the nine eigenvalues are close to each other, whereas for N > N^, 
three out of the nine eigenvalues become much larger than the others. 


( |3.3| ) although we do not introduce a cutoff in the temporal direction such as (^). In Fig. I 
(Left) we plot the results against N. At small N, it is almost independent of N. However, 
for N > Nc = 112, it begins to increase linearly with N. Figure |I| (Right) shows the N 
dependence of the expectation values (Aj (t)) of the nine eigenvalues of Tjj (t) evaluated at 
t = tpeak, where R'^{t) becomes maximum.^ For small N, there is no significant gap between 
the nine eigenvalues, whereas for N > Nc, we observe a big gap between (A 3 (tpeak)) and 
(A 4 (tpeak))- We will see in section |5| that the SO(9) symmetry is broken down to SO(3) 
after a critical time similarly to the original Lorentzian type IIB matrix model. 


4. Properties of the bosonic model for N < Nc 


In this section we discuss the properties of the bosonic model for N < Nc- In order to 


extract the time-evolution, we need to determine the block size n to be used in eq. (2.15). 
In Fig. ^ (Left) we plot the magnitude of the off-diagonal elements of Ai against the 
time separation Oq — for N = 110. The origin in the horizontal axis corresponds to the 
diagonal elements. We observe a nice scaling behavior for all the matrix elements. However, 
the magnitude falls off rather smoothly as one goes in the off-diagonal direction, which 
means that the dominant matrix configurations do not have a band-diagonal structure. 


In this situation, we cannot naturally define the block matrices (2.15) representing 
the state at each time and hence the notion of time-evolution becomes obscure. Let us 
nevertheless try to extract the “time-evolution” using n = 14 as the block size, which is 
the value obtained for N = W = 112 in the way described in the next section. In Fig. ^ 
(Right) we plot the expectation values (A* (t)) for N = 110. It turns out that there is only 
little t-dependence, and there is no clear gap between the eigenvalues for all t. 


'*In order to define lx{t) and Tij {t), we have to specify the block size n to be used in eq. ( [2.15[ ). See 
sections ^ and for the actual values of n used to obtain the results in Fig. 0 (Right). 
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Figure 2: (Left) The magnitude K^Oabl^ off-diagonal elements of Ai is plotted against 

the time separation aa — ctb for N = 110. (Right) The expectation values (A^ (t)) of the nine 
eigenvalues of Tij (t) are plotted against t for = 110. The block size is chosen to be n = 14. 



Figure 3: The extent of space (t) is plotted against t for iV = 64, 96 and 110. The block size 
is chosen to be n = 14 for all N. 

The situation for smaller N is similar to the N = 110 case. In Fig. |3|we plot the extent 
of space R‘^{t) as a function of t for N = 64, 96 and 110 obtained with the same block size 
n = 14. The dependence on N turns out to be modest. 

5. Properties of the bosonic model for N > 

In this section we study the properties of the bosonic model for N > N^. In Fig. ^ (Left) 
we plot the magnitude of the off-diagonal elements of Ai for N = 128. We find that 
the magnitude decreases rapidly as one goes away from diagonal elements. Moreover, the 
magnitude scales only for sufficiently large \aa — ctb\- From this observation, we identify 
the block size n as the number of points in the region where the off-diagonal elements do 
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Figure 4: (Left) The magnitude K^Oabl^ off-diagonal elements of Ai is plotted against 

the time separation —at, for N = 128. The scaling is observed only for sufficiently large \aa — at,|. 
(Right) The expectation values (At (t)) of the nine eigenvalues of Tij (t) are plotted against t for 
N = 128 with the block size n = 20. 

not scale. (In the present N = 128 case, we obtain n = 20. See below for more detail.) 
Using the block size n determined in this way, we can obtain the time-evolution. In Fig. ^ 
(Right) we plot the expectation values (A* (t)) for N = 128. In contrast to the situation for 
N < Nc, we observe the spontaneous symmetry breaking from SO(9) to SO(3) at a critical 
time tc similarly to the original Lorentzian type IIB matrix model.® 

In order to study the large-scaling property, we perform simulation for N = 256, 
384, 512 as well. In Fig. |5| (Top-Left) we zoom up the region near the origin in Fig. ^ 
(Left). From this figure, we determine the block size for N = 128 to be n = 20. Similarly, 
from the other figures in Fig. |^, we determine the block size for N = 256, N = 384 and 
512 to be n = 24, n = 28 and 32, respectively. 

The definition of the critical time tc is ambiguous at finite N. See Fig. ^ (Left), 
where we plot the expectation values (A* (t)) of the eigenvalues of Tij (t) against t for 
N = 512. Note first that the appearance of a gap between {X^it)) and {Xi{t)) signals the 
spontaneous symmetry breaking of SO(9) to SO(3). Let us therefore define the separation 
~ {Xj+i{t)). Then we find that the symmetric phase can be characterized by 
di{t) > d 2 {t) > ■ ■ ■ > d%{t), while in the broken phase we find d 2 {t) < d 3 {t). Therefore we 
may define the critical time tc by the largest value of t' such that di{t) > d 2 {t) > ■ ■ ■ > ds{t) 
holds for t < t'. For instance, the critical time tc obtained in this way for N = 512 from 
Fig. P (Left) is tc = —0.76559(7). Similarly we obtain tc = —0.76930(7) for N = 384. 
Applying the same procedure to smaller N, we find that the large-A^ scaling behavior in 
Fig. ^ (Right) becomes less clear due to finite N effects. We absorb these finite N effects 
by adjusting the value of tc slightly® from the one determined from the above procedure. 

®The fact that the spatial dimensionality after the spontaneous symmetry breaking turned out to be the 
same as in the original model is understandable from the view point of the mechanism suggested in ref. [44], 
which involves only the boson part of the action. 

®For N = 256, we shift by two data points and use tc = —0.82166(6) instead of tc = —0.76987(6). 
Similarly, for N = 128, we shift by four data points and use tc = —0.89472(7) instead of tc = —0.75798(7). 
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Figure 5: (Top-Left) The zoom up of the region near the origin in Fig. ^ (Left). We find 20 points 
in the region in which the scaling behavior is violated. Analogous plots for N = 256, N = 384, 
N = 512 are shown in the other panels, where we find 24, 28, 32 points in the non-scaling region, 
respectively. 


As is proposed in the original Lorentzian type IIB matrix model [44], we use the extent of 
space R{tc) at the critical time to fix the scale. Explicit values of R{tc) are given in table 
1^ together with the block size n and the critical time tc for each N. 

In Fig. ^ (Right) the extent of space R‘^{t) is plotted against t. The large-scaling 
behavior is observed by shifting the time coordinate so that the critical time comes to the 
origin and by plotting dimensionful quantities in units of R{tc). The observed large- 
scaling shows that the theory one obtains in the large-A1 limit is characterized by one scale 
parameter R{tc) and it does not contain any dimensionless parameters. 

It turns out that the behavior of R?{t) at t > tc can be fitted to an exponential 
function only for a finite range. At later times, it can be fitted well by a straight line, 
which corresponds to the power-law expansion 

R {t) oc . (5.1) 

Note that this behavior is observed within the scaling region, which implies that the sug¬ 
gested power law persists in the large-A^ limit at least for some time region. In Appendix]^ 
we present the results for the (5-|-l)D version of the bosonic type IIB matrix model. While 
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Figure 6: (Left) The expectation values (Ai (t)) of the nine eigenvalues of (t) are plotted against 
t for iV = 512 with the block size n = 32. (Right) The extent of space {t) normalized by (tc) is 
plotted against x = {t — tc) /R (tc) for N = 128, 256, 384 and 512. See table || for the values of the 
block size n, the critical time tc and the extent of space R{tc) at the critical time, which are used to 
make this plot. The solid line is a fit of the N = 512 data to R^ (t) jR? {tc) = a+ (1 — a) exp {bx) for 
1.0 < X < 1.85, which gives a = 0.9957(5) and b = 4.03(7). The dashed line is a fit of the N = 512 
data to R^ (t) jR? (tc) = cx + d for 1.85 <x< 2.5, which gives c = 34.3(6) and d = —55(1). 


we observe qualitatively the same behaviors, there are also some interesting quantitative 
differences. 

In order to understand the observed large-scaling further, we investigate how the 
continuum limit and the infinite volume limit in the temporal direction are achieved in the 
large-A^ limit. Here we restrict ourselves to > 256 since N = 128 is too close to the 
critical value Nc = 112. As the “lattice spacing” in the temporal direction, we consider 
the separation of data points in Fig. ^ (Right) in the horizontal direction. This quantity 
is actually t-dependent, and it can be defined more explicitly as ^ , where 6t is the 
difference of (1^ between adjacent blocks. In Fig. (Left) we plot this t-dependent 
“lattice spacing”, choosing the horizontal axis to be the same as in Fig. ^ (Right). We find 
clear tendency that the “lattice spacing” at the same point on the horizontal axis decreases 
as N increases. As the “volume” in the temporal direction, we define 


A = 


(^peak a) 

R{tc) 


(5.2) 


Using this quantity, we can also dehne an average “lattice spacing” e = A/u, where u is 
the number of data points within the region [tc, tpeak]- The values of e and A obtained 
for each N are given in table 2. We find that the average “lattice spacing” £ decreases and 
the “volume” increases as N becomes large. In Fig. ^ (Right) we plot e and A against 
N in the log scale. The straight lines represent fits to the power-law behaviors, although 
the behaviors may be subject to a qualitative change at larger N. In particular, it is an 
interesting dynamical question whether A —>■ oo or A —> const, in the large-A limit. In 
the former case, the expansion of space continues forever, since the time tpeak a-t the peak 
cannot be reached within a hnite time. In the latter case, on the other hand, the space 


- II - 














0.12 




Figure 7: (Left) The “lattice spacing” is plotted against {t — tc)/R{tc)■ (Right) The average 
“lattice spacing” e and the “volume” in the temporal direction A is plotted against N in the log 
scale. The straight lines represent fits to the power-law behaviors e = aN~P, where a = 0.20(1), 
p = 0.16(1) and A = bN'^, where b = 1.0(2), q = 0.18(3). 


N 

n 

tc 

R{tc) 

e 

A 

128 

20 

-0.89472(7) 

0.39270(2) 

— 

— 

256 

24 

-0.82166(6) 

0.30045(3) 

0.08297(2) 

2.7380(7) 

384 

28 

-0.76930(7) 

0.26580(3) 

0.07823(2) 

2.8943(6) 

512 

32 

-0.76559(7) 

0.24578(3) 

0.07417(2) 

3.1150(7) 


Table 1: The block size n, the critical time tc and the extent of space R{tc) at the critical time, 
which are used to make the plot in Fig. || (Right). We also present the explicit values of the average 
“lattice spacing” e and the “volume” A in the temporal direction, which are plotted in Fig. ^ 
(Right). 


stops expanding in a finite time and starts to shrink. By addressing this issue in the original 
model, one can, in principle, predict the fate of our Universe. 

6. Summary and discussions 

In this paper we have studied a bosonic model, which can be obtained from the Lorentzian 
type IIB matrix model by omitting the fermionic matrices. Due to the attractive potential 
between the eigenvalues of Aq arising from integrating out the spatial matrices at one 
loop, the eigenvalue distribution of has a finite extent even if one does not introduce 
the temporal cutoff (p.7|). In the original model, this attractive potential is canceled by 
the repulsive potential arising from integrating out the fermionic matrices at one loop and 
from the van der Monde determinant. (The simplified model for early times studied in 
ref. [57] has the same feature since the Pfaffian is replaced by the one-loop contribution.) 
Therefore, one would naively think that supersymmetry is playing an important role in the 
properties of the model that enable the extraction of a sensible time-evolution. 

Indeed for N < Nc, we observe that the extent of the eigenvalue distribution of 
is almost independent of N, and the dominant matrix configurations do not allow the 
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extraction of a sensible time-evolution. However, for N > N^, we find that the extent of 
the eigenvalue distribution of Aq grows linearly in N, and a sensible time-evolution can be 
extracted. We find that the SO(9) symmetry is spontaneously broken down to SO(3) at 
some critical time similarly to the original model [44]. The quantity such as R{t) shows a 
clear large-A^ scaling behavior. The growth of R{t) oc at late times is reminiscent of 
the behavior of the Friedmann-Robertson-Walker universe in the radiation dominated era. 

Let us recall that in the simplified model for early times, the growth of R{t) was 
observed to be exponential [57]. In that model, only the first term in ( 2.13| ) was used 
to represent the effect of fermionic matrices. We consider that the exponential expansion 
occurs also in the original model at early times as is suggested by direct Monte Carlo 
studies up to < 24 [58]. At late times, however, the subleading term in ( 2.13| ) becomes 
important due to the expansion of space, and that would affect the expanding behavior. 
Note that the repulsive potential for the eigenvalues of Aq is obtained from integrating out 
the fermionic matrices without the subleading term. Therefore one of the effects of the 
subleading term would be to make the repulsive potential less effective. Considering that 
the bosonic model mimics such a situation, we speculate that the exponential expansion 
in the original model changes into a power-law expansion at some point in time, where the 
subleading term in ( 2.131 ) becomes important. According to this scenario, the number of 
e-foldings is determined dynamically in the Lorentzian type IIB matrix model. It would be 
interesting to confirm the transition directly by simulating the original model. An attempt 
in doing this with a systematic approximation is in progress. 
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A. Details of Monte Carlo simnlation 


In this section we give the details on how we perform Monte Carlo simulation of the bosonic 
model (| 


First the delta functions in (3.3) are replaced by Gaussian potentials as 


1 


f"pot = X7C i^Tr (F;,.F^") + -7 l wrTr (A)" - 1 


N' 


1 


1 


N' 


(A.l) 


where the coefficients jc and 'Jl are taken large enough to fix each observable to the 
specified value. The values used in actual simulation are given in table 
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Another important issue we have to take care of is the spontaneous breaking of the 
shift symmetry Aq i—>• Aq + a\. For instance, let us consider calculating the expectation 
value B? (t) defined in (2.17). The peak of this quantity measured for each configuration 
fluctuates considerably. This reflects the ambiguity in choosing the origin of the time 
coordinate, and we should fix it before taking the ensemble average. Here we fix it by 
introducing a potential 


ksym — 


1 


-7syn 




- left 


— [Tr(Ai2 
N I ^ ' 


- right 


Tr {Ai) 

Tr(AT^ 


- left 


- right 


E E . 

i=l a-\-b<N+l 

E E . 

i=l a-\-b>N+l 


(A.2) 

(A.3) 

(A.4) 


where the values of the coefficient 7 sym used in our simulation are given in table | 
To summarize, the model we simulate is given by 


N 


z = 


daa dAi 


3 —‘S'efF 


a=l 


2 = 1 


•S'eff = -2 log A (a) + Tpot + Hsym . 


(A.5) 


The simulation of the model (A..5) can be performed by using the Hybrid Monte Carlo 
(HMC) method. First we rewrite the model by introducing auxiliary variables pa and 
{Xi)ab{a, 6 = 1,..., N) with the action 




(A.6) 


Here pa are real variables, whereas Xi are traceless Hermitian matrices. We update all 


the variables in the model (A. 6 ) in the following way. First we regard pa as the conjugate 


momenta of «« and Xi as the conjugate momenta of A,. Then we regard 5 hmc ia (|A. 6 D as 
the Hamiltonian H and solve the classical equations of motion obtained as the Hamilton 
equations 


dota 

II 

^1 CD 

II 

s 

dpa _ 

dH 


dT 

dT 

daa 

daa 

dAi 

dH 

- - X*, 

dXi 

dH 


dT 

dXi 

dT 

dA, 

dAi 


(A.7) 


for some fictitious time r. This part of the algorithm is called the Molecular Dynamics. 


In order to solve the Hamilton equations (A.7) numerically, we discretize them using the 
so-called leap-frog discretization, which maintains reversibility with respect to r. Starting 
from the previous conhguration at r = 0 , we obtain a new configuration at r = Tf by 
solving (AA) with the step size At so that Tf = W • At, where Nj- is the number of 
steps. We accept the new configuration with the probability min ( 1 , exp (— AS’hmc)); where 
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^‘S'hmc = 5'hmc ("^f) ~ 5'hmc (0), following the idea of the Metropolis algorithm to satisfy 
the detailed balance. The important point here is that 5'hmc is nothing but the Hamiltonian 
H, which is preserved in the classical dynamics if the equations (|A.7D are solved exactly. 
In fact, A5hmc becomes non-zero due to the discretization, but it is guaranteed to be a 
small quantity of the order of (Ar)^. By choosing sufficiently small At, the acceptance 
rate can be kept reasonably high, which enables the system to move around efficiently 
in the configuration space. Note also that the auxiliary variables pa and appear 

only as the Gaussian terms in ( A. 6 ). Therefore, we can update them independently by 
using normalized Gaussian random numbers. This procedure of refreshing the conjugate 
momenta should be done each time we start a Molecular Dynamics procedure. 

To summarize, the HMC algorithm as applied to our system can be described as follows. 


1. Generate initial configurations of Pa{0) and Aj (0) with Gaussian distribution oc 

e~2^a(Pa) and ^ respectively. 

2. Evolve the fields pa (r), Xi (r), aa (r) and Ai (r) for fictitious time Tf according to 
the discretized Molecular Dynamics. 

3. Accept the configuration of (rf) and A, (rf) obtained at the end of Molecular 

Dynamics with the probability min where AH = H (rf) — H {0). 


In the HMC algorithm, there are two parameters^ Ar and Tf. In the present work we 
choose them as in table 


N 

7c /iV^ 


Tsym 

Nr 

Ar 

trajectories 

128 

1 

100 

200 

20 

0.0015 

2 , 000,000 

256 

1 

100 

200 

10 

0.0008 

1,600,000 

384 

1 

100 

2,000 

10 

0.0004 

1 , 000,000 

512 

1 

100 

6,000 

10 

0.00025 

2,250,000 


Table 2: The values of the parameters yc, 7 l and ^sym in (A.5) used in our simulation. We also 
give the values of the parameters in the HMC algorithm: the number of steps Nr in the Molecular 
Dynamics and its step size Ar. In the last column, we give the number of “trajectories”, which 
represents how many times we solve the Molecular Dynamics after thermalization to achieve the 
statistics of our data. 


B. Results for the (5-|-l)D version of the bosonic model 

In this section we present our results® for a bosonic model that can be obtained by omitting 
fermionic matrices in the (5-|-l)D version of the type HB matrix model. The latter model 
is obtained formally by dimensional reducing the 6 D M = 1 super Yang-Mills theory to a 

^These parameters can be optimized as follows. For fixed rf, it is optimal to choose Ar so that Ar x 
(acceptance rate) is maximized. Then rf can be optimized to minimize the auto-correlation time in units 
of one step in the Molecular Dynamics. 

^Preliminary results shown in Fig. ^ (Right) are published in the proceedings [58]. 
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N 

7c/iV^ 


Tsym 

Nr 

Ar 

trajectories 

64 

1 

100 

1,000 

10 

0.0015 

200,000 

96 

1 

100 

1,000 

10 

0.001 

400,000 

128 

1 

100 

2,000 

10 

0.001 

2,400,000 


Table 3: The values of the parameters 7 c, 7 l and 7sym in (A.l) used in our simulation of the 
(5+l)D model. We also give the values of parameters in the HMC algorithm, (See caption of table 
0 for explanation.) 


point, and it consists of six bosonic matrices (^ = 1 ,..., 6 ) and four fermionic matrices 
'I'a (« = 1) • • • )4) representing four components of a 6 D Weyl spinor. The form of the 
bosonic part of the action is the same as that of the original type IIB matrix model, which 
is given in (p. 2 |). 



Figure 8: (Left) The extent (;^Tr (Aq)^) of the eigenvalue distribution of Aq is plotted against 
N for the (5+l)D model. It starts to increase a.t N = N^. = 34. (Right) The expectation values 
{Xi (t)) of the five eigenvalues of T^- (t) at t = tpeak are plotted against N for the (5+l)D model. 


In Fig. ^ (Left), we plot the extent (^Tr(Ao)^) of the eigenvalue distribution of Aq 
against N. In Fig. ^ (Right), we plot the expectation values A* (t) of the five eigenvalues 
of Tij (t) at t = tpeak against N. While the qualitative behaviors are the same as in the 
(9+l)D case shown in Fig. 2, we find that the critical N^. is smaller and the slope of the 
linearly increasing behavior of (;^Tr (Aq)^) for N > Nc is larger. We can understand this 
difference by considering the attractive potential between the eigenvalues of ^0 discussed 
below eq. ( p.l2 ). For general spatial dimensionality d, one obtains a factor (a) from 
integrating out the spatial matrices Ai at one loop. This factor acts as an attractive 
potential between the eigenvalues of Aq, and it is stronger for larger d. 

Below we discuss the properties of the (5+l)D model for N > Nc- (The parameters 
used in the simulation are listed in table |3|.) We have determined the block size to be 
n = 8 , 10, 12 for N = 64, 96, 128, respectively, from the fall-off of the off-diagonal elements 
of Ai as is done for the (9-|-l)D case in section In Fig. ^ (Left) we plot the expectation 
values {Xi (t)) of the five eigenvalues of Tij (t) for N = 128, which shows that the SO(5) 
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N 

n 

tc 

R{tc) 

e 

A 

64 

8 

-0.7248(5) 

0.1575(4) 

0.2281(6) 

1.825(5) 

96 

10 

-0.7692(3) 

0.1276(3) 

0.2157(4) 

2.157(4) 

128 

12 

-0.8037(1) 

0.1070(1) 

0.2048(2) 

2.457(2) 


Table 4: The block size n, the critical time tc, the extent of space R{tc) at the critical time, which 
are used in the (5+l)D model to make the plot in Fig. || (Right). We also present the explicit 
values of the average “lattice spacing” e and the “volume” A in the temporal direction, which are 
plotted in Fig. ^ (Right). 



t 



Figure 9: (Left) The expectation values (Ai (t)) of the five eigenvalues of Ty (t) are plotted against 
t for the (5+l)D model with N = 128. The critical time is determined as tc = —0.8037 (1). 
(Right) The extent of space {t) is plotted against x = {t — tc) /R{tc) for N = 64,96 and 128 
in the (5+l)D model. The solid line represents a fit of the N = 128 data to R^ (t) jR? {R) = 
a + (1 — a) exp (bx) for 0.4 < x < 1.2, which gives a = 0.839(9) and b = 2.91(6). The dashed line 
represents a fit of the N = 128 data to R^ (t) jR? {tc) = cx + d for 1.2 < x < 2.0, which gives 
c = 15.6(5) and d = —13.0(8). 


symmetry is broken spontaneously down to SO(3) after a critical time. From this kind of 
figures, we can determine the critical time tc for each N as described® in section |^. 

In Fig. 1^ (Right) we show the large-scaling behavior of the extent of space {t). 
Explicit values of R{tc), together with the block size n and the critical time tc, which are 
used to make this plot, are given in table Q. The power-law expansion (^) is observed at 
late times similarly to the (9-|-l)D model. 


In Fig. 10 (Left) we plot the t-dependent “lattice spacing”, which shows how the 
continuum limit is achieved as N increases. The average “lattice spacing” e and the 
“volume” A in the temporal direction are given in table In Fig. (Right) we plot them 
in the log scale. The straight lines represent fits to the power-law behaviors. 


^Unlike in the (9+l)D case, there was no need to adjust the value of tc to obtain the large-A scaling 
behavior in Fig. ^ (Right). 
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Figure 10: (Left) The “lattice spacing” is plotted against {t — tc)/R{tc) for the (5+l)D 

model. (Right) The average “lattice spacing” £ and the “volume” in the temporal direction A is 
plotted against N in the log scale for the (5+l)D model. The straight lines represent fits to the 
power-law behaviors £ = aN~P, where a = 0.44(2), p = 0.16(1) and A = 6 A®, where b = 0.30(2), 
q = 0.43(1) using all the data. 
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